version 16

clear
clear mata
clear matrix
set more off
set matsize 11000
set maxvar 30000
cap log off
capture log close
pause on

cap ssc install outreg2
cap ssc install reghdfe
cap ssc install ftools


global path "/Users/abehrer/Desktop/"

global code "$path/Replication code/Code/"
global data "$path/Replication code/Data/"
global output "$path/Replication code/Output/"
global Table "$output"

	use "$data/Replication_data.dta", replace
	
	local a_dry = " prec_dry below40_dry FND_Max40to45_dry FND_Max45to50_dry FND_Max50to55_dry FND_Max55to60_dry FND_Max60to65_dry FND_Max75to80_dry  FND_Max65to70_dry "
	
	local a_wet = " prec below40_wet FND_Max40to45_wet FND_Max45to50_wet FND_Max50to55_wet FND_Max55to60_wet FND_Max60to65_wet FND_Max75to80_wet  FND_Max65to70_wet "
	
	local b_dry = " c.Dividends#c.prec prec    c.Dividends#c.below40_dry c.Dividends#c.FND_Max40to45_dry c.Dividends#c.FND_Max45to50_dry c.Dividends#c.FND_Max50to55_dry c.Dividends#c.FND_Max55to60_dry c.Dividends#c.FND_Max60to65_dry c.Dividends#c.FND_Max75to80_dry  below40_dry FND_Max40to45_dry FND_Max45to50_dry FND_Max50to55_dry FND_Max55to60_dry FND_Max60to65_dry FND_Max75to80_dry FND_Max65to70_dry c.Dividends#c.FND_Max65to70_dry"	
	
	
	local b_wet = " c.Dividends#c.prec prec    c.Dividends#c.below40_wet c.Dividends#c.FND_Max40to45_wet c.Dividends#c.FND_Max45to50_wet c.Dividends#c.FND_Max50to55_wet c.Dividends#c.FND_Max55to60_wet c.Dividends#c.FND_Max60to65_wet c.Dividends#c.FND_Max75to80_wet  below40_wet FND_Max40to45_wet FND_Max45to50_wet FND_Max50to55_wet FND_Max55to60_wet FND_Max60to65_wet FND_Max75to80_wet  FND_Max65to70_wet c.Dividends#c.FND_Max65to70_wet"	
	
	
	cap mat define tempbinfigure=J(20,6,.)
	local row=1
	
	cap mat define tempdivfigure=J(48,6,.)
	local row2=1
	
	//Dry bulb regression
	rename (days90plus_dry FND_Max85to90_dry FND_Max80to85_dry  ) (days90plus FND_Max85to90 FND_Max80to85 )

	
	reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85  `a_dry', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
			est store m1
			local raw=_b[days90plus]
			foreach var in below40_dry FND_Max40to45_dry FND_Max45to50_dry FND_Max50to55_dry FND_Max55to60_dry FND_Max60to65_dry FND_Max65to70_dry FND_Max70to75_dry FND_Max75to80_dry FND_Max80to85 FND_Max85to90 days90plus{
				cap mat def tempbinfigure[`row',1]=_b[`var']
				cap mat def tempbinfigure[`row',2]=_se[`var']
				cap mat def tempbinfigure[`row',3]="`var'"
			local row=`row'+1
			} 
			g rawdayimpact90=`raw'
	
	rename (days90plus FND_Max85to90 FND_Max80to85 ) (days90plus_dry FND_Max85to90_dry FND_Max80to85_dry  ) 	

			
			//Create Variables for future projections using dry bulb temps
		rename (days90plus_dry FND_Max85to90_dry FND_Max80to85_dry  ) (days90plus FND_Max85to90 FND_Max80to85 )

			rename Dividends x
				foreach n in pc_Div pc_DivL1 mean_Div4 mean_Div2 { 
					rename `n' Dividends
					label var lny_nonag "`n'"
					label var Dividends "Dividends"
					
					reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85 Dividends c.Dividends#c.days90plus c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85 `b_dry', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
					
					est store `n'dry	
						g reg`n'_dry_B1=_b[days90plus]
							local direct`n'_dry=_b[days90plus]
						g reg`n'_dry_B2=_b[c.Dividends#c.days90plus]
							local indirect`n'_dry=_b[c.Dividends#c.days90plus]
						g reg`n'_dry_effect=reg`n'_dry_B1+(c.Dividends*reg`n'_dry_B2)
					
					
					rename Dividends `n'
				}
				
		//Regressions for figure 1
			foreach n in mean_Div4  { 
					rename `n' Dividends
					label var lny_nonag "`n'"
					label var Dividends "Dividends"
					
					reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85 Dividends c.Dividends#c.days90plus c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85 `b_dry', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
					foreach var in FND_Max80to85 FND_Max85to90 days90plus c.Dividends#c.FND_Max80to85 c.Dividends#c.FND_Max85to90 c.Dividends#c.days90plus {
				cap mat def tempdivfigure[`row2',1]=_b[`var']
				cap mat def tempdivfigure[`row2',2]=_se[`var']
				cap mat def tempdivfigure[`row2',3]=e(df_r)
				local row2=`row2'+1
			} 
					
					rename Dividends `n'
				}
				
				
		rename (days90plus FND_Max85to90 FND_Max80to85 ) (days90plus_dry FND_Max85to90_dry FND_Max80to85_dry  ) 
	g dayimpact90=`directmean_Div4_dry'
	g interactimpact=`indirectmean_Div4_dry'
	
		rename x Dividends	
	
	//Wet bulb regression	
	rename (days90plus_wet FND_Max85to90_wet FND_Max80to85_wet  )  (days90plus FND_Max85to90 FND_Max80to85 ) 

	
	
	reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85   `a_wet', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
			est store m2
			local raw=_b[days90plus]
			local row=1
			foreach var in below40_wet FND_Max40to45_wet FND_Max45to50_wet FND_Max50to55_wet FND_Max55to60_wet FND_Max60to65_wet FND_Max65to70_wet FND_Max70to75_wet FND_Max75to80_wet FND_Max80to85 FND_Max85to90 days90plus{
				cap mat def tempbinfigure[`row',4]=_b[`var']
				cap mat def tempbinfigure[`row',5]=_se[`var']
				cap mat def tempbinfigure[`row',6]="`var'"
			local row=`row'+1
			} 
			
		//Temperature bin figure for dry and wet bulb temperatures
		preserve	
			clear
			
			svmat2 tempbinfigure

	
		g tbin=_n
		drop if tbin>12
		
			replace tempbinfigure1=tempbinfigure1*100
			replace tempbinfigure2=tempbinfigure2*100
			replace tempbinfigure4=tempbinfigure4*100
			replace tempbinfigure5=tempbinfigure5*100
	
			g HI95_dry=(1.96*tempbinfigure2)+tempbinfigure1
			g LI95_dry=tempbinfigure1-(1.96*tempbinfigure2)

			g HI99_dry=(2.06*tempbinfigure2)+tempbinfigure1
			g LI99_dry=tempbinfigure1-(2.06*tempbinfigure2)

			g HI95_wet=(1.96*tempbinfigure5)+tempbinfigure4
			g LI95_wet=tempbinfigure4-(1.96*tempbinfigure5)

			g HI99_wet=(2.06*tempbinfigure5)+tempbinfigure4
			g LI99_wet=tempbinfigure4-(2.06*tempbinfigure5)
			
			mvencode H* te* L*, mv(0)


			twoway (rcap HI99_dry LI99_dry tbin, lc(gs10) ) ///
					(line tempbinfigure1 tbin, lc(ebblue) lp(solid) ) ///
					(rcap HI99_wet LI99_wet tbin, lc(gs10) lp(dash)  ) ///
					(line tempbinfigure4 tbin, lc(edkblue) lp(longdash) ///
					xti("Temperature Bin") xlabel(1 "Below 40" 2 "40-45" 3 "45-50" 4 "50-55" 5 "55-60" 6 "60-65" 7 "65-70" 8 "70-75" 9 "75-80" 10 "80-85" 11 "85-90" 12 "Above 90", angle(45)) ytitle("Change in annnual payroll (%)") legend(off) yline(0))
					
					graph export "$output/tempbin_figure.png", as(png) replace
		restore			
			
	
	

		//Regressions with various measures of Dividends
		local i=3	
		rename Dividends x
		foreach n in pc_Div pc_DivL1 mean_Div4 mean_Div2 { 

			rename `n' Dividends
			label var lny_nonag "`n'"
			label var Dividends "Dividends"
			
			reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85  c.Dividends#c.days90plus c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85 Dividends`b_wet', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
			
			est store m`i'
			local i=`i'+1
			rename Dividends `n'
		}
		
		//Regressions for figure 1
		local i=3	
		foreach n in mean_Div4 { 

			rename `n' Dividends
			label var lny_nonag "`n'"
			label var Dividends "Dividends"
			
			reghdfe lny_nonag days90plus FND_Max85to90 FND_Max80to85  c.Dividends#c.days90plus c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85 Dividends`b_wet', absorb(i.fips i.year i.state#c.year) cluster(i.fips) 
			local row2=1
			foreach var in FND_Max80to85 FND_Max85to90 days90plus c.Dividends#c.FND_Max80to85 c.Dividends#c.FND_Max85to90 c.Dividends#c.days90plus {
				cap mat def tempdivfigure[`row2',4]=_b[`var']
				cap mat def tempdivfigure[`row2',5]=_se[`var']
				cap mat def tempdivfigure[`row2',6]=e(df_r)
				local row2=`row2'+1
			} 
			
			
			rename Dividends `n'
		}
		
	
	rename (days90plus FND_Max85to90 FND_Max80to85 x) (days90plus_wet FND_Max85to90_wet FND_Max80to85_wet  Dividends) 
	label var Dividends "Dividends"
	
	g days90plus=0
	g FND_Max80to85=0
	g FND_Max85to90=0
	label var FND_Max85to90 "Days \$85-90\degree\$F"
	label var FND_Max80to85 "Days \$80-85\degree\$F"
	label var days90plus "Days \$>90\degree\$F"
	
	
	//Figure 1
		preserve	
			clear
			
			cap graph drop graph1 
			cap graph drop graph2
			
			svmat2 tempdivfigure
			
			
	
		g tdiv=_n
		drop if tdiv>12
			replace tdiv=1.3 if tdiv==2 
			replace tdiv=1.6 if tdiv==3
			replace tdiv=4.3 if tdiv==5 
			replace tdiv=4.6 if tdiv==6
		
		g tdiv2=tdiv+.025
		replace tdiv=tdiv-0.025
	
			g HI95_dry=(1.96*tempdivfigure2)+tempdivfigure1
			g LI95_dry=tempdivfigure1-(1.96*tempdivfigure2)

			g HI99_dry=(2.06*tempdivfigure2)+tempdivfigure1
			g LI99_dry=tempdivfigure1-(2.06*tempdivfigure2)

			g HI95_wet=(1.96*tempdivfigure5)+tempdivfigure4
			g LI95_wet=tempdivfigure4-(1.96*tempdivfigure5)

			g HI99_wet=(2.06*tempdivfigure5)+tempdivfigure4
			g LI99_wet=tempdivfigure4-(2.06*tempdivfigure5)
			
			g tstat_dry=abs(tempdivfigure1/tempdivfigure2)
			g tstat_wet=abs(tempdivfigure4/tempdivfigure5)
			
			g pval_dry=2*ttail(tempdivfigure3,abs(tstat_dry))
			g pval_wet=2*ttail(tempdivfigure6,abs(tstat_wet))
			
			g tstat_dry1=tstat_dry if pval_dry>0.05
			g tstat_wet1=tstat_wet if pval_wet>0.05
			
			g tstat_dry2=tstat_dry if pval_dry<=0.05 & pval_dry>0.01
			g tstat_wet2=tstat_wet if pval_wet<=0.05 & pval_wet>0.01
			
			g tstat_dry3=tstat_dry if pval_dry<=0.05 & pval_dry>0.01
			g tstat_wet3=tstat_wet if pval_wet<=0.05 & pval_wet>0.01
			
			g tstat_dry4=tstat_dry if pval_dry<=0.01 & pval_dry>0.001
			g tstat_wet4=tstat_wet if pval_wet<=0.01 & pval_wet>0.001
			
			g tstat_dry5=tstat_dry if pval_dry<=0.001 & pval_dry>0
			g tstat_wet5=tstat_wet if pval_wet<=0.001 & pval_wet>0
			
			
			
			mvencode H* te* L* t* p*, mv(0)


			twoway (rcap HI99_dry LI99_dry tdiv in 1/3, lc(gs10)  ) ///
					(scatter tempdivfigure1 tdiv in 1/3, mc(ebblue) ms(S)  ) ///
					(rcap HI99_wet LI99_wet tdiv2 in 1/3, lc(gs10) lp(dash)   ) ///
					(scatter tempdivfigure4 tdiv2 in 1/3, mc(edkblue) ms(D) ) ///
					(bar tstat_dry1 tdiv in 1/3, yaxis(2) barwidth(0.05) color(edkblue) lc(none) ) ///
					(bar tstat_wet1 tdiv2 in 1/3,yaxis(2) barwidth(0.05) color(edkblue) lc(none) ) ///
					(bar tstat_dry3 tdiv in 1/3, yaxis(2) barwidth(0.05) color(eltblue) lc(none) ) ///
					(bar tstat_wet3 tdiv2 in 1/3,yaxis(2) barwidth(0.05) color(eltblue) lc(none) ) ///
					(bar tstat_dry4 tdiv in 1/3, yaxis(2) barwidth(0.05) color("245 184 208") lc(none) ) ///
					(bar tstat_wet4 tdiv2 in 1/3,yaxis(2) barwidth(0.05) color("245 184 208") lc(none) ) ///
					(bar tstat_dry5 tdiv in 1/3, yaxis(2) barwidth(0.05) color("232 114 161") lc(none) ) ///
					(bar tstat_wet5 tdiv2 in 1/3,yaxis(2) barwidth(0.05) color("232 114 161") lc(none)  ///
					xti("Temperature Bin", height(6)) xlabel( 1 "80-85" 1.3 "85-90" 1.6 "Above 90") ytitle("T-stat			Change in annnual payroll (%)", margin(0 -3 0 0)) ytitle("",axis(2)) legend(off) yline(0) ysc(noline) ylabel(-0.0020 "0" -0.00183 "2" -0.00166 "4" -0.001 "-0.1" -0.0005 "-0.05" 0 "0" 0.0005 "0.05", nogrid) yscale(range(-0.002 0.0005) axis(1)) ylabel(none, axis(2)) yscale(range(0 30) axis(2)) name(graph1))
					
					graph export "$output/tempdiv_figureA.png", as(png) replace
					
				twoway (rcap HI99_wet LI99_wet tdiv in 4/6, lc(gs10) ) ///
					(scatter tempdivfigure4 tdiv in 4/6, mc(edkblue) ms(D) ) ///
					(rcap HI99_dry LI99_dry tdiv2 in 4/6, lc(gs10) lp(dash)  ) ///
					(scatter tempdivfigure1 tdiv2 in 4/6, mc(ebblue) ms(S)) ///
					(bar tstat_wet1 tdiv in 4/6, yaxis(2) barwidth(0.05) color(edkblue) lc(none)) ///
					(bar tstat_dry1 tdiv2 in 4/6,yaxis(2) barwidth(0.05) color(edkblue) lc(none)) ///
					(bar tstat_wet3 tdiv in 4/6, yaxis(2) barwidth(0.05) color(eltblue) lc(none)) ///
					(bar tstat_dry3 tdiv2 in 4/6,yaxis(2) barwidth(0.05) color(eltblue) lc(none)) ///
					(bar tstat_wet4 tdiv in 4/6, yaxis(2) barwidth(0.05) color("245 184 208") lc(none)) ///
					(bar tstat_dry4 tdiv2 in 4/6,yaxis(2) barwidth(0.05) color("245 184 208") lc(none)) ///
					(bar tstat_wet5 tdiv in 4/6, yaxis(2) barwidth(0.05) color("232 114 161") lc(none)) ///
					(bar tstat_dry5 tdiv2 in 4/6,yaxis(2) barwidth(0.05) color("232 114 161") lc(none) ///
					xti("Dividends by Temperature Bin") xlabel( 4.6 `""Dividends x" "80-85""' 4.3 `""Dividends x" "85-90""' 4 `""Dividends x" "Above 90""')  ytitle("T-stat		Change in annnual payroll (%)", margin(0 -1 0 0)) xsc(reverse) ytitle("",axis(2)) legend(pos(6) col(2) row(2) holes(1 2 3 6 7) order(2 "Dry Bulb" 4 "Heat Index" - "t-statistic:" 5 "p>0.05" 7 "p<0.05" 9 "p<0.01" 11 "p<0.0001")) yline(0) ysc(noline) ylabel(-0.0005 "0" -0.00038 "2" -0.00025 "4" 0 "0" 0.001 "0.1" 0.0005 "0.05" 0.0015 "0.15", nogrid) yscale(range(-0.0005 0.0015) axis(1)) ylabel(none, axis(2)) yscale(range(0 30) axis(2)) name(graph2))
					
					graph export "$output/tempdiv_figureB.png", as(png) replace
			grc1leg graph1 graph2, legend(graph2)
					
					graph export "$output/tempdiv_figure.png", as(png) replace
					
		restore
		
	
	
	//Main table
	esttab m1 m2 mean_Div4dry m5 using "$Table/main_table.tex", replace f ///
		label booktabs b(5) p(5) eqlabels(none) alignment(S S) mtitles("\shortstack{Dry-Bulb}" "\shortstack{Heat Index}" "\shortstack{Dry Bulb//5-year avg.}" "\shortstack{Heat Index//5-year avg.}" )  ///
		keep(days90plus  FND_Max85to90 FND_Max80to85 Dividends c.Dividends#c.days90plus  c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85) ///
		star(* 0.10 ** 0.05 *** 0.01) ///
		refcat(, nolabel) ///
		stats(N r2, fmt(%9.0fc 2) layout("\multicolumn{1}{c}{@}" "\multicolumn{1}{S}{@}") labels(`"Observations"' `"\(R^{2}\)"' `"County FE"' `"Month FE"' `"Year FE"' `"State $\times\$ Year"' `"State $\times\$ Year$^{2}\$"' `"State $\times\$ Year$^{3}\$"'))

	//SI Table with alternative specifications
	esttab pc_Divdry m3 pc_DivL1dry m4 mean_Div2dry m6 using "$Table/main_table_SI.tex", replace f ///
		label booktabs b(5) p(5) eqlabels(none) alignment(S S) mtitles("\shortstack{Dry Bulb//Current}"  "\shortstack{Heat Index\\Current}" "\shortstack{Dry Bulb\\Lagged}" "\shortstack{Heat Index\\Lagged}"  "\shortstack{Dry Bulb\\11-year avg.}" "\shortstack{Heat Index//11-year avg.}")  ///
		keep(days90plus  FND_Max85to90 FND_Max80to85 Dividends c.Dividends#c.days90plus  c.Dividends#c.FND_Max85to90 c.Dividends#c.FND_Max80to85) ///
		star(* 0.10 ** 0.05 *** 0.01) ///
		refcat(, nolabel) ///
		stats(N r2, fmt(%9.0fc 2) layout("\multicolumn{1}{c}{@}" "\multicolumn{1}{S}{@}") labels(`"Observations"' `"\(R^{2}\)"' `"County FE"' `"Month FE"' `"Year FE"' `"State $\times\$ Year"' `"State $\times\$ Year$^{2}\$"' `"State $\times\$ Year$^{3}\$"'))
	
	drop days90plus
	drop FND_Max80to85
	drop FND_Max85to90
	
	run "$code/Table2 Fig 2A 2B.do"
		
		
	rename regmean_Div4_dry_effect regmean_Div4_effect 
	
		preserve
			collapse (mean) regmean_Div4_effect rawdayimpact90, by(PovDecile)	
			save "$data/ForProjections_Admin1.dta", replace
		restore
		
		preserve
			collapse (mean) regmean_Div4_effect rawdayimpact90, by(PovCentile)	
			save "$data/ForProjections_Admin2.dta", replace
		restore	
		
		run "$code/Table 3.do"
		
		erase "$data/ForProjections_Admin2.dta"
		erase "$data/ForProjections_Admin1.dta"
	
	run "$code/Fig 1A 1B.do"
	
	run "$code/SI Table 1 Pre2003.do"
	run "$code/SI Table 2 Top code.do"
	run "$code/SI Table 3 AC.do"
	run "$code/SI Table 4 HE.do"
	run "$code/SI Table 6 Dry HI Comp.do"
	//Summary stats, SI figures
	run "$code/SI Sum Stats.do"
		
	clear	

		
		
